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NOMENCLATURE . 


A, A' 

matrix defined in equations (16) and (17) 
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Planck's function (eq.(2b)) 
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quantity defined in equations (15b) and (A12) 
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speed of light 


E 
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exponential integral of order n 


F 

nondimensional radiative flux (eq. (6a)) 


f- 

boundary function 
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integral of the kernel function 


9 

boundary function defined in. equation (7b) 


H 

integral of B over n (eqs. (15a) and (A15)) 


h 

Planck's constant, or function defined in equation (7b) 


I 

one-half the total step number 


J 

extremal function 


K 

kernel function 


k 

Boltzmann's constant 


L 

physical length 


M 

number of iterations 


N 

number of spectral modes 
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radiative heat 

* 

S 

source function 

• 

s 

s 

nondimensional source function, — ~ 
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sgn 

r , 1+ if X > 0 

sgn(x> = ^ <0 


T 

temperature 


X 

dummy variable 

iii 
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k 
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X 
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4 

B 

J 

J 

k 


independent variable, altitude 
variable, 

absorption coefficient defined in equation (3a) 
delta function, (eq. (6a)), or variational of the function 
increment of step function 
increment of frequency, 

f.y 

independent distance variable, * J 6 dz/ 
optical thickness for the kth frequency, a^n 

T 

normalized temperature, 

1 ^ 

one-half the thickness, J* 3 dy 
mass absorption coefficient 
constant associated with integral equation 
frequency, or an arbitrary function (eq. (A7)) 

nondimens ional frequency, 

mass density 

Stefan-Boltzmann constant 

optical thickness defined in equation (4b) 

variational function 

solution 

Subscripts 

at the lower boundary 
at the upper boundary 
at Jth point 

indexes of summation, or at ith or jth step point 
frequency interval corresponding to kth mode 


IV 
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St the boundary 
frequency 

even or odd function 


V 

+ 


+ 


positive or negative 


Superscript 

direction 
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AN ITERATIVE SOLUTION OF AN INTEGRAL EQUATION FOR RADIATIVE TRANSFER 

BY USING VARIATIONAL TECHNIQUE 


Kenneth K. Yoshikawa 
Ames Research Center 


SUMMARY 


An effective iterative technique is introduced to solve a nonlinear inte- 
gral equation frequently associated with radiative transfer problems. The 
problem is formulated in such a way that each step of an iterative sequence 
requires the solution of a linear integral equation. The advantage of a pre- 
viously introduced variational technique which utilizes a stepwise constant 
trial function is exploited in this report to cope with the nonlinear problem. 
The method is simple and straightforward. Rapid convergence is obtained by 
employing a linear interpolation of the iterative solutions. Using absorption 
coefficients of the Milne-Eddington type, = ayg(p,T), which are applicable 
to some planetary atmospheric radiation problems, solutions are found in terms 
of temperature and radiative flux. These solutions are presented numerically 
and show excellent agreement with other numerical solutions. 


INTRODUCTION 


The variational approach is attractive for solving the integral equations 
of kinetic theory bound ary -value problems and aero-gasdynamics . It is versa- 
tile, accurate, fast, and particularly useful in solving linear integral equa- . 
tions (refs. 1, 2, and 3). The variational function may be accurately 
represented by polynomials which lead to extremely difficult and complex 
evaluations. To alleviate these difficulties, reference 2 examined a simple 
step-function representation. For nonlinear integral equations (e.g., 
Fredholm's type) the number of solution techniques available is extremely 
limited and numerical calculation times are long (ref. 4) . 

In this report a new iterative technique is introduced to solve a nonlin- 
ear integral equation. The problem is formulated in such a way that each step 
of an iterative sequence requires the solution of a linear integral equation. 

In a previous paper, an efficient variational technique (using a stepwise con- 
stant trial function) was described for the solution of the second kind of 
linear Fredholm's equation. Consequently, the advantages of the variational 
method for linear equations can be exploited within the framework of an itera- 
tive sequence for the nonlinear problem. The results of the present method 
are compared with other numerical solutions and the significant advantages of 
the present method are demonstrated. 



ANALYSIS 


Radiative Transfer Equation 


The basic transfer equation for one-dimensional atmospheric radiation is 
(e.g. , refs. 3 and 4) 

00 _ CO 

= .2^ j E2(T^)dv - 2pJ q-^ E2(t^^ - 


/vL 


27t J* PK^ dv f S^(2 ')£'i(|t^ - t\)dt 


+ 47T 


rj pK^B^CDdv 

0 


The net radiative heat flux q is the direct integral of equation (1), 

/ “ oo 

- 2Trj* - x^)dv 

fi •'n 


where 


i«” (> \>L 

IttJ dv j B^(iZ’)sg'n(T^ - t)E2(jT^ - t|)dt 

O' *0 


,y 


f (y) ~ I PK d = T (L) 

V J vw vB 
0 


^„(T) =/ x-^ 


dx 


B^(T) = (2;zv3/c2)/(e'^''/^ - 1) 


( 1 ) 


(2a) 


(2b) 


The symbols q* and q~ are the boundary fluxes due to the surface emissivity 

(e) and reflectivity (r) at boundaries A and B, respectively (see sketch (a)). 
The symbol ic^ is the spectral absorption coefficient assumed to be of the 
Milne-Eddington type; that is. 


P<y = * B(P,T) (3a) 

where the strength is approximated by the "picket fence" model referring to 
N different rectangles, 

“v "" “fe '^k - ^k+l ^ 2,...N) (3b) 
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The use of the Milne-Eddington type 
model, however, maybe limited to the 
planetary atmospheres where pressure- 
induced transition is important, as 
for example, in the atmospheres of 
Saturn, Uranus, and Neptune. The geom- 
etry and coordinate system under con- 
sideration are shown in sketch (a) . 

To simplify the analysis, symmetric 
coordinates have been used, a pro- 
cess which leads to the following 
transformation for optical thickness 



q (y) 


Boundary B 


where 


T = T, 

V k 




0 




(4a) 


(4b) 




y = 0 TTTTT'TTTTTTTTTT 
1 


Boundary A 


Sketch (a) . - Geometry and coordinate 
system. 


The heating rate in equation (1), expressed in terms of the nondimensional 
internal source function e(p,2!), is assumed to be given by 


d<7 = 


Hp- 


on 


(4c) 


so that equation (1), when combined with equations (3) and (4), can be rewrit- 
ten as 

N j- 
k=i L 




( 5 ) 


where s, 6, n are the nondimensional internal source function, (assumed to be 
known), ratio of atmospheric temperature to the boundary temperature T/T^ (to 
be determined), and optical thickness, respectively, and 
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B^(6) = ttJ 6(a^)B^(e)dv/(a!r^) = e'^JCv/G) 

0 

00 

^Ak ' 

00 

•'o 

00 

J(v/0) = 6Ca^)-(v/e)3d(v/e)/^e‘''/®- 1^ 


V = hv/kr 


B 


6(a^) 


= (o 


for the frequency 


otherwise 


Nondimensional radiative flux F is given by 

N 


F(n) i «/aI* =Vp^(,p 
fe=l 


where 


> (6a) 


(6b) 


r^k . 

+ 2J sgniT]j^ - n')B;;,(e)£' 2 (|n^ - n'|)dn' (6c) 

~^k 

Note that the black body radiation has the relation of (which includes S, for 

= 0 ) 

N 

2 B^(9) = s'* (6d) 

k=i 


where 6 = T/T^. 

Equation (5) is a nonlinear integral equation to be solved for 9 = 9(n). 
The description of the basic steps in solving the nonlinear equation by the 
present method may be briefly outlined as follows: First, the nonlinear inte- 

gral equation is rearranged to take the standard form of Fredholm's integral 
equation (second kind). Second, with a pair of arbitrarily selected tempera- 
ture distributions (stepwise-constant, initially isothermal distributions for 
the present problem) the nonhomogeneous part of the integral equation is eval- 
uated analytically for the pair of temperature distributions. The two result- 
ing linear integral equations are then solved by the variational method 
introduced previously, using a stepwise -const ant trial function. Finally, an 
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iterative method for calculating a new pair of temperature distributions is 
described; this method uses the linear interpolation (or extrapolation) between 
the pair of previous temperature distributions and its resulting linear solu- 
tions. The iterative procedure is continued until sufficient convergence is 
attained. 


Linear Representation 

For the present problem, equation (5) is recast into the standard form of 
Fredholm's integral equation (second kind) as 


/ ®'l 

B^(n’)£i(|ni - n'|)dn' 

-£l 


(7a) 


where 

N 

/(n^) = j (5) + ai^i(b^) ^ 

k=2 

i I i 

\ = 


(7b) 


Note that the / function in equations (7b) can be evaluated analytically since 
temperature is represented by a stepwise constant function. Initial guesses 
are two arbitrary isothermal approximations, for example, boundary temperature 
6=8. and 6 = O-. The detailed iteration procedure is discussed in Iterative 

Calculations and Discussions (also in appendix B) . The reference function 
a segment of B given by Eq. (6a), can be selected arbitrarily from the other 
spectral functions B^. However, it is desirable to pick Bj from the larger 
value of calculated for the isothermal conditions, so that more accurate 

solutions can be maintained. 

The term containing the source function will be omitted from the present 
calculation. However, the general procedure will be exactly the same as for 
the rest of the terms in the f function. 
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Extremal Function 


It can be shown (e.g., ref. 1) that the linear integral equation 
(Fredholm's equation of the second kind) containing the symmetric singular ker- 
nel function X(n,n'), that is, 

'f'Cn) = /(n) + k X(n,n')'J^(n')dn' (8) 

can be derived from the functional 

= if-si ^(n,n')^Cn')dn' - 2/(n)Jdn (9) 

where J takes a minimum value at (}) = ij;, which is the solution of equation (8). 


For the simplicity of numerical calculations, the linear trial function c|), 
the boundary function /, and the extremal functional tJ((t)) are divided into 
even and odd functions, designated by subscripts + and - respectively, as 


where 


(j> = + (p_ 

/ = 4 ^ 

4 = J [/(n) + /(-n)] 

4 = J [fM - /C-n)] 


(10a) 

(10b) 

(10c) 


(lOd) 


Note that the extremal function due to the cross product of an even and odd 
function has been eliminated since it is zero. Thus, the variational of the 
function J becomes 


= 0 


( 11 ) 


The number of points to be solved is significantly reduced (by one-half) 
by equation (11) since the even and odd function of cj) can be determined inde- 
pendently. Following equation (9), the extremal function J associated with 
equations (7) can be written as 


J(<).) = J 


J* <t>^(nj)dnj ' <f'(nj)dTijJ* <|)(nj)Si(|nj -n'|)dn' 

~^i "^1 


- 2 


/■ 


4>(ni)9'^(n^)dni 





( 12 ) 
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where is replaced by <|). For the conditions given by equation (10) and 

(11), even and odd functions of <i> can be determined for the minimum value of 

J. 


Solution by Step Function 

The even and odd functions of a trial function <)> in equation (10a) are 
replaced by stepwise constant functions (see sketch (b)) given by 



Sketch (b).- Notations for the odd step 
function 


nCn) 

■i 


for ^ n < 


(13a) 


<l) (n) = 0 


for nj - n < 


where the interval is equally divided by 21 points. An = 2£/(2J - 1), and by 
definition 
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(13c) 


21 I 

^(A+)^- = 0 ; 2^ = 0 ; or (A_)_^ 

J=l J=i 


J-1 


- E 


A major advantage of using the A function (leather than the direct 0. function) 
is that many terms pertaining to the integral of the kernel function can be 
eliminated from the calculations by virtue of equation (13c) (see appendix A). 
Even though the system defined in sketch (b) is for only the odd function, simi- 
lar notations apply for the even function, with the exception that the func- 
tion is symmetric and the value at the center (<j>^) is not necessarily zero. 
After some algebraic manipulation, substitution of equations (13) into equa- 
tion (12) yields the following equations (see appendix A). 


where 


C^(^ + ) = 


I 

E 

i=l 


r . T 2 


E fv’, 


J=i 




I I 

i=i j=i 

k=i 


(14) 


N - I 


z=i 

21 


- lE IE'M^ 

k=2 ' i=i j=i 


* 2 

i=i ^ 






^%-h 

= {.G . 

1,3 

--i*3k 



= [-^al 

'■\h - 







= (An)^ 

E 


/ — V 
V 


- ^^3-i^k ’ = 

^^2ih - 0 


(i - J) I 


(15a) 


/ 
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(15b) 


{B .), = B, (0 .) 


and 






*\ “ '’Ic^l ■ \ 


Note that Bj^ in equation (15b) includes all frequency intervals for 
In equation (14), minimizing J with respect to A. yields the following 

For the even function <|) : 


where 


9J(6 ) 

3 = 1 


^ A. = 0 (i = 1,2,...J) 
J ^ 




A. . - A.. 

^3 J'Z' 


for 3 < i) 


A, = (4 


N 

3 1 * Z)Wi’ 


k=2 




^^.ih " ■ l^^Ak * ^Bk^^^^li^k ■ ^'^l,2J+l-i^fe^ 
21 

J=i 

- 2[iH.)^ - I 


(15c) 

= a^. 
results : 

(16a) 

(16b) 

(16c) 

(16d) 


(16e) 
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For the odd function <p_ : the index of summation is reduced to J - 1 by sub- 

stituting the relation given by equation (13c) into equation (14) 


I-l 


dJ(<p 
9 (A 


= 0 (i = 1,2,...J-1) 

--'i J = 1 ^0 3 ^ 


where 




= 4[(n^), - (n^.),] - 

* 7- forjii 


/I ' /I ’ 

1 . . = A .. 

-iQ QT' 


for j < i 


N 


and 


= W-), - 

(^^)l =-(^^1 - * ^^l,2J+l-f^l " " *^S,I+1^J 

/ 

21 

- (C, T 1 ) 7 ] + Z-7(AS.)7 [(G. O 7 + (G. , O 7 - (G^T-07 

^ IjJ+l-^fe-' J fc*- ^3 k ^ ^,2I■*■l-Q■’k '• I^-’k 


(17a) 


(17b) 


(17c) 

(17d) 


- * ^'‘2in-i\ - <-hh - 

Equations (16a) and (17a) and I and (J - 1) simultaneous linear equations for 
A., where coefficients of these equations are calculated from a given tempera- 
ture distribution. Successive approximations follow thereafter until the con- 
vergence in (J>^ or 6^ is obtained (see appendix B) . From equations (10a), 
(13a), (16a), and (17a) 







( 18 ) 
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Evaluation of the extremal function J given by equation (14) and the 
procedure for solving the essential equations (16) and (17) , are systematic 
and straightforward, requiring minimal computing time (Computing time is 
roughly proportional to I^NM, where N is total number of frequency intervals 
and M is a number of iterations.). Numerical calculations for J as low as 3 
show surprisingly good accuracy, and the computing time on an IBM 360/67 for a 
general case using J = 13 is nominal [i.e., less than 1/10 sec per frequency 
interval per iteration (see appendix B)]. ^ 

The radiative fluxes, from equations (6b) and (6c), are calculated by 

N N 

F(-£) = 

k=i k=i 

where 

21 

‘ ^Ak - 

J = 1 

21 

and 

I 

= 2(.F^ . [- Cn^O,] (19c) 

J = 1 

However, as noted in reference 2, the radiative flux F computed by the present 
method is not exactly constant. The mean value, suggested in reference 2, is 
therefore employed here; 

F = J [F(-£) + f(0)] (20) 


ITERATIVE CALCULATIONS AND DISCUSSIONS 


For nonlinear equations, the successive substitution or successive 
approximation methods are extremely time consuming (ref. 4) and frequently 
experience stability problems if the initial approximation is not adequately 
selected. The efficiency of the calculation depends on simplifications made 
in the analysis, the desired order of accuracy, convergence rate, and the 
scope of the results of interest (i.e., whether interest lies in the macro- 
scopic physical quantity or solution itself). The convergence rate also 
depends on the function to be solved; for example, in the present problem, the 
solution in terms of temperature 0 converges much faster than the direct 
solution of the function ij) itself. For a faster convergence, the successive 
approximation method used for the present problem is therefore modified by 
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considering a linearly weighted effect of two independent initial guesses and 
their corresponding solutions. That is, two sets of arbitrarily selected 
temperature distributions, 0i and 62 , are assumed (Initially, isothermal dis- 
tributions are used.). The integrated Planck's function, Bi (0) and the 
solution function <() are then computed from equations (15b) and (18), respec- 
tively. The results are then used to obtain one new temperature distribution 
03 by linear interpolation (or extrapolation) and a second new temperature 
distribution ©4 by linear interpolation between ©3 and either 0 i or 02 (see 
appendix B) . With these new data, the iteration procedure is repeated until 
the temperatures 0 converge. The present results show that the fourth or 
fifth iterate provides sufficient convergence. 

Sample calculations for several models of absorption coefficients have 
been selected for this study. Included in these is the model from reference 4 
which is used to compare the present results with other numerical solutions. 

It is noted the number of spectrals N (where B > 1) in the absorption coef- 
ficient has no significant effect on the convergence of the solution in the 
present method. However, the selection of the reference Planck function 
(e.g., if much smaller is selected) will obviously affect the accuracy of 
the solution. 


The results for temperature distribution and normalized flux using the 
simple spectral model from reference 4 are presented in tables 1 and 2. In 
these tables the optical thickness is selected such that Tjt = 1 and the num- 
ber of steps is varied from 3 to 13 in order to determine the effect of step 
size on the computed radiative flux. The present results show surprisingly 
good agreement with the numerical values of reference 4 even for the minimum 
number of steps, 1=5. (The midpoint values of each step function are tab- 
ulated in the illustrations.) Table 1 shows the results of reference 4 and 
the present results for 1=3 and 13, so that comparisons can be made for 
coincident optical distance n. The results of normalized radiative flux (F) 
are presented in table 2 for various step sizes. The step size for J = 13 is 
more than sufficient to retain good accuracy for the cases considered in this 
paper. 

The present method has been used to calculate several other selected 
sample cases using absorption coefficients described by a larger number of 
spectral modes and a variety of optical thicknesses. The absorption coeffi- 
cient and temperature 0 for a spectral model consisting of 36 nonzero and 
nonrepeated values (model 3) are presented in tables 3 and 4. The optical 
thickness selected for table 4 is = 2. Note that the calculation will be 
significantly reduced when the absorption coefficient contains zero or 
repeated values. The results for a three step mode (I = 3) shows good agree- 
ment with those of larger steps (J = 13) for temperature distribution and 
radiative flux (Note: The last decimal in table 4 may be not correct even 

though the solution converges to the fifth decimal place, since the computer 
carries only five significant figures.). Results of radiative flux using 
model 3 with several optical thicknesses are presented in table 5 (J = 13) . 

An interesting feature of the variational method is that accurate values 
for a physical quantity, such as radiative flux F for the present case, can 
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be obtained without having an accurate solution for 0. For example, a linear 
(or in some cases, only a one step constant mode) trial function provides 
remarkably accurate values for physical quantities even though the solutions 
themselves are not necessarily accurate (ref. 1) . The present result substan- 
tiates this conclusion; in other words, the results for macroscopic physical 
quantities show much faster convergence than the solution itself. Conse- 
quently, a large number of steps and/or iterations are not usually necessary. 
Examples of this convergence on the solution (temperature at n = -0.96 and 
n = 0.96) and radiative flux are illustrated in table 6. 

The boundary gas temperature, even though it is not considered in this 
paper, may be calculated in two ways: (1) by using a small step interval at 
the boundaries; and (2) evaluating the right-hand side of equation (7a) by 
substituting the present step solution for at n = ±H (also at any point) . 

An iterative method used in conjunction with the variational technique 
which uses a step-function provides the following advantages: 

1. Accurate descriptions of the solutions can be obtained with a simple 
iteration scheme (nonlinear case) . 

2. Rapid convergence of the physical quantities may be obtained without 
accurate solutions. 

3. The procedure of numerical calculations is straightforward. 

4. Nonlinear problems which are difficult to handle analytically by 
other methods (e.g., the polynomial trial function or successive approximation 
methods) may be exploited by a method similar to the present technique. 

CONCLUDING REMARKS 

With the aid of the variational method to the linear integral equation, 
an iterative technique has been applied to solve the nonlinear integral equa- 
tion which is frequently involved in physical and gas dynamic problems. It 
has been demonstrated that the iterative technique, in conjunction with the 
variational method which uses the step function as a trial function, provides 
effective and accurate descriptions of the solutions to these problems and 
related average physical quantities. 

The method is relatively simple, straightforward, and requires nominal 
computing time. It is expected that the method can be further extended for 
solving other nonlinear integral equations of radiative transfer or problems 
in kinetic theory. 

Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., January 22, 1973 
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APPENDIX A 


EVALUATION AND SIMPLIFICATION OF INTEGRATIONS 
INVOLVING THE KERNEL FUNCTION 


Integrations of kernel functions using a step function are usually 
relatively simple. Further simplifications which are important in reducing 
the computation process are realized by eliminating all terms from the inte- 
grated function when these terms vanish in certain summations. 


Integrals of Symmetric Kernel Functions 


Some major results involving the integrals of the symmetric functions can 
be summarized. For the exponential integrals 


'n 

Similarly 

Jl 


J ffiCln - n' |)dn = sgn(x\^ - n')[-l + £’2(l\ - p' hi + 1 - - n') 

•'n . 


®2(|n - n' - n')|^- y + Esdn^ - n’ |) ♦ j - SsC* - i') 


•'n • 


sg?2(n* - n-)dn' = “|n- " n-| + - n.) 

u "u J u 


sg?z(n' - n^)£’2(|n' - n^|)dn' = " n^-l) ~ - \) 


CM) 

(A2) 

CA3) 

(A4) 


Simplification of Integrated Functions 

Integration of the symmetric kernel function with a singularity required 
in equation (12) if is a step function can be evaluated by equation (Al) 
through (A4) as follows: 




n - n ' 1 )dn dn ' 


= (z - n^) + (z - n^O - |n^ - n^-| 



+ £’3(£, - nO + £^3 (z - n-) 
^ j 


(AS) 
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Since the symmetric matrix G.. of equation (AS) will be summed over i and/or j 

'I'O 

for determining the extremal values of J (eq. all terms in G.. except 

(n- - n_0 will vanish from the summation as a result of the conditions imposed 
for the step function by equation (13c). Thus, G can be redefined as 




(A 6 ) 


Integration of Symmetric Functions 


Since the extremal function of equations (9) or (12) is evaluated for 
both even and odd boundary or trial functions (eqs. (9)), some simplifications 
with regard to the single and double integrations may be made; assuming K 
represents the symmetric kernel, 

£ £ 

f <j)+(n)dn f v(n' ) [if(n,n' ) ± ^(-n,n')]dn' 


= f v(n')dn' f <J)+(n) [^(n,n') ± Z(-n,n’)]dn 

£ £ 

= 2 r <fi^(n)dn f v(n')-S^(ri,ri')dn (A 7 ) 

•'-0 - 


where v(n) is an arbitrary function. Also 


Z 

f 4>+(ri) [^(Z,T]') ± K(Z,-n)]dri - 2J i}>^ (n) [A:(£,n) ± ^(«.,-n)]dTi 1 (A8) 

•'-£ - -Z ~ ! 


Equations (A7) and (A8) have been incorporated in the results indicated by 
equation (14) . 


Reduction of the Index of the Summation 


Since the solution is divided into the even and odd functions, the index 
of double summation appearing in equation (14) can be rearranged using equa- 
tions (13) as follows: 

f E".- "i.- = E ‘iE 

i=i i=i j=i 


J I 

" S * ^2J+l-i,2J+l-j + ^i,2J+l-j + ^2J+l-i, p 

i=l j = l 


J 


= E 

i=i .1=1 


A . (G . . 




(A 9 ) 
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where 


G • • — G • • 

^3 J'Z' 


2 1*1-3 ~ ^2I+l-i,j 


CAIO:) 


Note that the minus sign in the last equation (A9) is for the even function 
(or A ) while the plus sign is for the odd function (A). The results of equa- 
tion t-A9) have been used in obtaining equation (14) of the text. 


From the condition given by equation (13c) , the index of summation for the 
odd function' in equation (A9) can be further reduced. (For simplicity, the 
minus sign indicating the odd function has been dropped from the following 
equation. ) 


I I 

EA. y^A .(G. . + G. , .) 

3 ^3 ^J2J+l-J-' 

i=l 3=1 


I 

I-l 



^ A . (G .. + G._,j.. . - G . — - 
^ 3 13 t,2I+l-j rl 


i=i 

3=1 


I-l 

J-1 



^ ' A.(G.* + G, -T-. • - Q _ 

^ 3^ ^3 21 + 1-3 ^il 

^i,I+l 

i=l 

3 = 1 


- G . j. 

3 I 

— + (7 7-7- + Gy. r 1 ) 

JJI+1 II JJI+1-' 

(All) 


2 J -*1 3 been replaced by G ^3 T^*l 

ion (All) is given by equation (17b) 


The symbols Gj.^ and 
The result of equation” (All) 


^J,J+1 CAIO)), 

given by equation (17b) of the text. 


Miscellaneous Integration 

Other integrated functions appearing in equation (14) and (15) are derived 
and denoted as follows: 



(n)E2(^ + n)dn = 


21 

^ (A+)^ f G2C*- + n)dn 


)^[^3C^ + - ff3(2il)] 


By equation (13c), 


21 


21 

E (AJ .C . 

i=l 


E (4.),- 


21 


^3(2^) = ff3(2il)^(A^)^ = 0 
i=l i=l 


(A12) 
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It 


that 


and 




H ^3 * n .) 


Since thf, 

t\"'at ""Side the 


/. . f 

z^I " j 
21 

^i(n) 

t^l z- 


^(■0;dr 


^Cejdr 




bracJcet i 


the 


»7 


iast 


fAjj) 


WI4; 




^.hUhaa 


CAIS) 
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APPENDIX B 


THE METHOD OF ITERATION 


Since equation (5) and (7) are nonlinear integral equations, a modified 
successive substitution method using a pair of initial temperature distribu- 
tions is introduced for solving the resulting equations; in other words, the 
boundary function f in equation (7) or matrices A and A' in equations (16) 
and (17) are first calculated from the initial set of arbitrarily selected 
temperature distribution 6(1) and solved for (j> at each n • • 

The value 4> obtained from equation (18) is matched with the value (B) i 
calculated by equation (6) for the given 6(1). The same procedure is repeated 
for the second set of arbitrarily selected temperature distribution 6(2). From 
this pair of initial temperature distributions, 6(1) and 6(2), and their 
resulting values, (B) j and (p, the first new temperature distribution 6(3) is 
linearly interpolated at each point. The second new temperature 6(4) is then 
calculated by the linear interpolation (or extrapolation) between 6(3) and 
either 6(1) or 6(2). The scheme of the linear interpolation is formulated as 
follows (for simplicity, the subscript 1 pertaining to the reference function 
has been omitted) : 


Temperature 6(3) is given by 

. . (6^(2) - eJ(D] [♦ (!). B.(D] 

IV2) - b!(2).;,(1)4[1)) 


Temperature 6(4), depending on the conditions, is given by 


(a) if 

[<1,^(3) - B^ (3)] [6^(3) - ej(2)] > 0 

4 4 [9^(3) - ej(i)][<^.(i) - B (1)] 

^ ~ ' [^,-(3) - ^,-(3) - <^,-(1) + B,-(D] 


or 


(Bl) 


,(B2) 

\ 



(b) (replacing all the indexes 1 in equation (B2) with the index 2), 


e;(4) 


6 *!( 2 ) 


[6^(3) - 6^(2) ] (2) - B^(2)] 

[<^^(3) - B^(3) - <p^(2) + B^(2)) 


(B3) 
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For the next iterations, 0(5) and 0(6), the indexes 1 and 2 in equations (Bl) 
through (B3) are replaced by the indexes 3 and 4. Then the successive itera- 
tion M follows until (|) converges to B, as 

e.iM) - B.iM - 1) ^ 

^rp) — - ^ 0T(W) ^ 

For the present problem, the fifth iteration (M = 5) provides sufficient accu- 
racy. Note that the two initial temperatures selected for the present compu- 
tations are the isothermal distributions 0(1) = 0^ and 0(2) - where 0^ and 
0^ represent the respective temperatures at each of the boundaries. 

From the results of the present calculations (by J = 13), the average 
computing time required for the solution and flux is less than 1/10 sec per 
spectral mode per iteration on an IBM 360/67; in other words, for N = 36 and 
M = 5, the computing time was about 18 sec. 
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TABLE 1.- TEMPERATURE DISTRIBUTION: 6. = 0.5 , 0„ = I.O e = = 1.0, 

A D Ad 


= 0, AND = 1.0 


Model 1: 

ai = 1.0(vi = 0) 
«2 = 0.5 (v2 = 3) 

Model 2 : 

ai = 0.5(vi = 0) 
02 = 1.0(V2 = 3) 




Present 

method 


Present 

method 


Reference 




Reference 



n 

4 




4 






J = 3 

I = 13 


1=3 

J = 13 


0.7408 






0.7095 

. 



-.4 

.7726 

0.7723 

0.7726 

.7432 

0.7428 

0.7432 

-.2 

■ .8148 


.8148 

.8148 

.7882 

.7881 

.7882 


.8484 


.8483 

.8484 

.8247 

.8245 

.8247 

.2 

.8781 


.8780 

.8781 

.8574 

.8573 

.8574 

.4 

.9066 , 


.9066 

.9065 

.8899 

.8899 

.8899 

.5 

.9235 


--- 

— 

.9098 

— 

— 


TABLE 2.- NORMALIZED FLUX (-F) ; BOUNDARY CONDITIONS SAME AS TABLE 1. 


Models 

Present method 

Reference 

4 

(X 

1 

a 

2 

I = 3 

1=8 

1 

J = 13 

1.0 

0.0 

0.7856 

0.7860 

0.7860 

0.7860 

1.0 

.1 

.7371 

— 

.7376 

.7376 

1.0 

.3 

.6637 

— 

.6644 

.6644 

1.0 

.5 

.6092 

— 

.6100 

.6101 

1.0 

1.0 

.5177 

.5187 

.5188 

.5188 


1.0 

.5351 

.5361 

.5361 

.5362 

.5 

1.0 

.5693 

.5701 

.5702 

.5702 

.1 

1.0 

— 

— 

.6427 

.6427 

.0 

1.0 

.6696 

.6703 

.6703 

.6703 
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TABLE 3.- ABSORPTION COEFFICIENT FOR MODEL 3. 


V 

Av 

a- 

V 

V 

Av 

a- 

V 

0 

0.25 

0.65 

4.05 

0.04 

0.82 

.25 

.75 

.10 

4.09 

.01 

.72 

1.00 

.50 

.80 

4.10 

.05 

.62 

1.50 

.25 

.25 

4.15 

.03 

.52 

1.75 

.05 

.35 

4.18 

.01 

.42 

1.80 

.02 

.88 

4.19 

.01 

.32 

1.82 

.03 

.78 

4.20 

.55 

.22 

1.85 

.05 

.68 

4.75 

.25 

.55 

1.90 

.10 

.58 

5.00 

.25 

.275 

2.00 

.10 

.48 

5.25 

.25 

.75 

2.10 

.15 

.38 

5.50 

.50 

.45 

2.25 

.25 

.60 

6.00 

.25 

.90 

2.50 

.50 

.20 

6.25 

.25 

.975 

3.00 

.25 

1.00 

6.50 

.25 

.70 

3.25 

.25 

.30 

6.75 

.50 

.15 

3.50a 

.50^^ 

.85^ 

7.25 

.25 

.95 

4.00 

.03 

.40 

7.50 

.25 

.50 

4.03 

.02 

.92 

7.75 

00 

.05 


^Reference frequency and absorption coefficient selected for k = 1. 
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TABLE 4.- TEMPERATURE DISTRIBUTION (6) AND NORMALIZED FLUX (-F) ; BOUNDARY 

CONDITIONS SAME AS IN TABLE 1. MODEL 3iN = 36), = 2.0. 

h 


Optical Thickness, 

Temperature, 

0 

1=5 

J = 13 

-0.96 


0.73340 

-.88 


.74897 

-.80 

0.76202 

.76184 

-.72 


.77330 

-.64 


.78380 

-.56 


.79358 

-.48 


.80281 

-.40 

.81174 

.81159 

-.32 


.81997 

-.24 


.82804 

-.16 


.83583 

-.08 


.84338 

.00 

.85075 

.85075 

.08 


.85792 

.16 


.86495 

.24 


.87185 

.32 


.87867 

.40 

.88530 

.88541 

.48 


.89209 

.56 


.89878 

.64 


.90548 

.72 


.91228 

.80 

.91915 

.91925 

.88 


.92657 

.96 


.93476 

Flux, -F 

.5706 

.57124 
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TABLE 5.- NORMALIZED FLUX (-F) AS A FUNCTION OF OPTICAL THICKNESS (x^): 

L/ 

BOUNDARY CONDITIONS SAME AS IN TABLE 1. MODEL 3: (N = 36, J = 13). 


Optical thickness. 

Normalized flux, 
-F 

0 

1.00000 

0.5 

0.79393 

1.0 

.69628 

2.0 

.57124 

3.0 

.49076 

4.0 

.43307 

5.0 

.38907 

6.0 

.35412 

7.0 

.32552 

8.0 

.30161 


TABLE 6.- CONVERGENCE OF THE SOLUTION (TEMPERATURE AND FLUXES): 
BOUNDARY CONDITIONS SAME AS TABLE 4. 


Number of iteration, 
M 

Temperature, 

0(n) 

Normalized flux, 
-F 

n = -0.96 

n = 0.96 

1 

0.5000 

0.5000 

- 

2 

1.0000 

1.0000 

— 

3 

.7493 

.9298 

0.5703 

4 

.7354 

.9339 

.5711 

5 

.7334 

.9348 

.5712 

6 

.7334 

.9348 

.5712 

7 

.7334 

.9348 

.5712 

8 

.7334 

.9348 

.5712 

9 

.7334 

.9348 

.5712 
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